HerbPt1 16S Figures

require(phyloseq)
require(tidyverse)
require(reshape2)
require(dplyr)
require(ggplot2)
require(vegan)

Load data order, factors, and create a mode (chemical, hand, non-treated) column.

ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])


index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])



ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])

ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])

create alpha diversity tables

alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md) 

Joining, by = "Barcode_ID_G"

alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Handweeded", "Atrazine-Mesotrione", "Dicamba", "Glyphosate"))

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df  Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   66537   16634   0.099  0.982
Residuals   51 8553352  167713               

aov_t2<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df   Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   652901  163225   0.779  0.545
Residuals   49 10272922  209651               

aov_t3<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df   Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   433435  108359   0.465  0.761
Residuals   50 11641251  232825               

remove outliers and rare reads with less than 2 total reads

ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_rare) != "G166SG")

ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)

ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.05, ps_trans)

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

Run 0 stress 0.1117941 
Run 1 stress 0.1117967 
... Procrustes: rmse 0.00125091  max resid 0.01384803 
Run 2 stress 0.1134329 
Run 3 stress 0.1117728 
... New best solution
... Procrustes: rmse 0.01276998  max resid 0.1385506 
Run 4 stress 0.1118649 
... Procrustes: rmse 0.01087397  max resid 0.1023859 
Run 5 stress 0.1117959 
... Procrustes: rmse 0.01364467  max resid 0.1486566 
Run 6 stress 0.1115948 
... New best solution
... Procrustes: rmse 0.004588589  max resid 0.05391228 
Run 7 stress 0.111766 
... Procrustes: rmse 0.00442201  max resid 0.05494826 
Run 8 stress 0.1117957 
... Procrustes: rmse 0.01142902  max resid 0.1397853 
Run 9 stress 0.1115974 
... Procrustes: rmse 0.0004616654  max resid 0.005151518 
... Similar to previous best
Run 10 stress 0.1172805 
Run 11 stress 0.1117979 
... Procrustes: rmse 0.01131905  max resid 0.1393504 
Run 12 stress 0.1117954 
... Procrustes: rmse 0.01201186  max resid 0.1419067 
Run 13 stress 0.1118117 
... Procrustes: rmse 0.01155646  max resid 0.1407441 
Run 14 stress 0.1117657 
... Procrustes: rmse 0.004425952  max resid 0.05523005 
Run 15 stress 0.1134278 
Run 16 stress 0.1117966 
... Procrustes: rmse 0.01230531  max resid 0.1482626 
Run 17 stress 0.1115972 
... Procrustes: rmse 0.0004124194  max resid 0.004505644 
... Similar to previous best
Run 18 stress 0.1115948 
... New best solution
... Procrustes: rmse 0.0001752304  max resid 0.00151237 
... Similar to previous best
Run 19 stress 0.1117956 
... Procrustes: rmse 0.01141319  max resid 0.1396299 
Run 20 stress 0.1117961 
... Procrustes: rmse 0.01230562  max resid 0.1477826 
*** Best solution repeated 1 times

ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1959743 
Run 1 stress 0.1958498 
... New best solution
... Procrustes: rmse 0.00736199  max resid 0.08432568 
Run 2 stress 0.1967565 
Run 3 stress 0.1969623 
Run 4 stress 0.195836 
... New best solution
... Procrustes: rmse 0.004183262  max resid 0.04410732 
Run 5 stress 0.1972166 
Run 6 stress 0.1972004 
Run 7 stress 0.1967507 
Run 8 stress 0.1972118 
Run 9 stress 0.1955227 
... New best solution
... Procrustes: rmse 0.01439727  max resid 0.1055675 
Run 10 stress 0.1967629 
Run 11 stress 0.1955951 
... Procrustes: rmse 0.01280822  max resid 0.1482606 
Run 12 stress 0.1970466 
Run 13 stress 0.1967622 
Run 14 stress 0.1955933 
... Procrustes: rmse 0.01278164  max resid 0.1478142 
Run 15 stress 0.1958344 
... Procrustes: rmse 0.01355004  max resid 0.1053769 
Run 16 stress 0.1960805 
Run 17 stress 0.1965852 
Run 18 stress 0.1955274 
... Procrustes: rmse 0.001218077  max resid 0.009525113 
... Similar to previous best
Run 19 stress 0.1955259 
... Procrustes: rmse 0.001157442  max resid 0.01064565 
Run 20 stress 0.1955251 
... Procrustes: rmse 0.0006025696  max resid 0.004705696 
... Similar to previous best
*** Best solution repeated 2 times

full_ord_rare<-ggordiplots::gg_ordiplot(ord = ord_rare, groups = data.frame(sample_data(ps_rare))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

full_ord_rare$plot + theme_classic()

ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", trymax= 300, maxit=1000)

Run 0 stress 0.07436383 
Run 1 stress 0.07356844 
... New best solution
... Procrustes: rmse 0.05257797  max resid 0.5365183 
Run 2 stress 0.07356423 
... New best solution
... Procrustes: rmse 0.05700177  max resid 0.5049506 
Run 3 stress 0.07475784 
Run 4 stress 0.08233187 
Run 5 stress 0.07371847 
... Procrustes: rmse 0.03580052  max resid 0.2670563 
Run 6 stress 0.0748452 
Run 7 stress 0.07432888 
Run 8 stress 0.07356475 
... Procrustes: rmse 0.00101334  max resid 0.01257088 
Run 9 stress 0.07351843 
... New best solution
... Procrustes: rmse 0.03301646  max resid 0.277146 
Run 10 stress 0.07481794 
Run 11 stress 0.08269937 
Run 12 stress 0.07412137 
Run 13 stress 0.07564989 
Run 14 stress 0.08308348 
Run 15 stress 0.07434076 
Run 16 stress 0.07435862 
Run 17 stress 0.07790567 
Run 18 stress 0.07413747 
Run 19 stress 0.07388433 
... Procrustes: rmse 0.05802289  max resid 0.4890175 
Run 20 stress 0.0721808 
... New best solution
... Procrustes: rmse 0.04835115  max resid 0.5098029 
Run 21 stress 0.07432931 
Run 22 stress 0.07439095 
Run 23 stress 0.0741388 
Run 24 stress 0.07215009 
... New best solution
... Procrustes: rmse 0.001801121  max resid 0.02249514 
Run 25 stress 0.07401767 
Run 26 stress 0.07207119 
... New best solution
... Procrustes: rmse 0.02931455  max resid 0.2673238 
Run 27 stress 0.07804942 
Run 28 stress 0.07293172 
Run 29 stress 0.07458269 
Run 30 stress 0.07630386 
Run 31 stress 0.08146302 
Run 32 stress 0.08008002 
Run 33 stress 0.07439049 
Run 34 stress 0.08424865 
Run 35 stress 0.07343365 
Run 36 stress 0.07356858 
Run 37 stress 0.07618202 
Run 38 stress 0.07475861 
Run 39 stress 0.07256438 
... Procrustes: rmse 0.02950194  max resid 0.267304 
Run 40 stress 0.08201316 
Run 41 stress 0.0866015 
Run 42 stress 0.07209557 
... Procrustes: rmse 0.001905361  max resid 0.02317822 
Run 43 stress 0.08455037 
Run 44 stress 0.08007292 
Run 45 stress 0.07433954 
Run 46 stress 0.07382298 
Run 47 stress 0.07444042 
Run 48 stress 0.07414279 
Run 49 stress 0.07353217 
Run 50 stress 0.0731359 
Run 51 stress 0.07398746 
Run 52 stress 0.07214984 
... Procrustes: rmse 0.02925501  max resid 0.2670558 
Run 53 stress 0.07357666 
Run 54 stress 0.4158152 
Run 55 stress 0.07395362 
Run 56 stress 0.08199542 
Run 57 stress 0.07565232 
Run 58 stress 0.07207392 
... Procrustes: rmse 0.001457162  max resid 0.01630803 
Run 59 stress 0.07352626 
Run 60 stress 0.07438306 
Run 61 stress 0.08008342 
Run 62 stress 0.07356382 
Run 63 stress 0.07388934 
Run 64 stress 0.07398744 
Run 65 stress 0.0735644 
Run 66 stress 0.07374194 
Run 67 stress 0.07372377 
Run 68 stress 0.07356866 
Run 69 stress 0.08315409 
Run 70 stress 0.07346563 
Run 71 stress 0.07432932 
Run 72 stress 0.07419749 
Run 73 stress 0.07371894 
Run 74 stress 0.07475716 
Run 75 stress 0.07356427 
Run 76 stress 0.07389694 
Run 77 stress 0.07356259 
Run 78 stress 0.07398762 
Run 79 stress 0.07357632 
Run 80 stress 0.07451023 
Run 81 stress 0.07358983 
Run 82 stress 0.0761558 
Run 83 stress 0.07353107 
Run 84 stress 0.07436762 
Run 85 stress 0.07421491 
Run 86 stress 0.07659576 
Run 87 stress 0.08662864 
Run 88 stress 0.07565735 
Run 89 stress 0.07293179 
Run 90 stress 0.07478225 
Run 91 stress 0.07371931 
Run 92 stress 0.07411583 
Run 93 stress 0.07381286 
Run 94 stress 0.07214946 
... Procrustes: rmse 0.02928908  max resid 0.2677357 
Run 95 stress 0.07214635 
... Procrustes: rmse 0.002841601  max resid 0.02883108 
Run 96 stress 0.08146367 
Run 97 stress 0.07209518 
... Procrustes: rmse 0.00213922  max resid 0.02322161 
Run 98 stress 0.07436743 
Run 99 stress 0.07436394 
Run 100 stress 0.07475812 
Run 101 stress 0.07353144 
Run 102 stress 0.07347476 
Run 103 stress 0.08005355 
Run 104 stress 0.07631888 
Run 105 stress 0.08008087 
Run 106 stress 0.07634158 
Run 107 stress 0.07390232 
Run 108 stress 0.08005734 
Run 109 stress 0.07388468 
Run 110 stress 0.07451177 
Run 111 stress 0.0722182 
... Procrustes: rmse 0.02933714  max resid 0.2670719 
Run 112 stress 0.07565095 
Run 113 stress 0.073713 
Run 114 stress 0.07727503 
Run 115 stress 0.07360437 
Run 116 stress 0.07207283 
... Procrustes: rmse 0.0002854145  max resid 0.003421485 
... Similar to previous best
*** Best solution repeated 1 times

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}

remove one sample with no vegetation measurement.

ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)

ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
ps_adonis(ps_dmn_57)

Ordination plots DMN

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1) 
T2_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")


ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1) 
T3_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")

Ordination plots rarefied

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491 
Run 1 stress 0.181468 
... Procrustes: rmse 0.05480179  max resid 0.2537838 
Run 2 stress 0.1812637 
... Procrustes: rmse 0.02578262  max resid 0.1092543 
Run 3 stress 0.1885435 
Run 4 stress 0.1823917 
Run 5 stress 0.1820307 
Run 6 stress 0.1873026 
Run 7 stress 0.1831937 
Run 8 stress 0.1877083 
Run 9 stress 0.1816042 
Run 10 stress 0.1821568 
Run 11 stress 0.1816494 
Run 12 stress 0.1815751 
Run 13 stress 0.1812635 
... Procrustes: rmse 0.02573893  max resid 0.109003 
Run 14 stress 0.180928 
... New best solution
... Procrustes: rmse 0.02695614  max resid 0.1109664 
Run 15 stress 0.1812552 
... Procrustes: rmse 0.01647177  max resid 0.08102514 
Run 16 stress 0.1814158 
... Procrustes: rmse 0.02898979  max resid 0.1057012 
Run 17 stress 0.1837645 
Run 18 stress 0.1833285 
Run 19 stress 0.1822759 
Run 20 stress 0.1812636 
... Procrustes: rmse 0.01524903  max resid 0.08041434 
Run 21 stress 0.1819126 
Run 22 stress 0.1837641 
Run 23 stress 0.1812579 
... Procrustes: rmse 0.01876556  max resid 0.09463233 
Run 24 stress 0.1826156 
Run 25 stress 0.1831936 
Run 26 stress 0.1820917 
Run 27 stress 0.181525 
Run 28 stress 0.180942 
... Procrustes: rmse 0.01975739  max resid 0.09890697 
Run 29 stress 0.1823225 
Run 30 stress 0.1854467 
Run 31 stress 0.1810715 
... Procrustes: rmse 0.03976112  max resid 0.1470678 
Run 32 stress 0.1916469 
Run 33 stress 0.1815252 
Run 34 stress 0.1816037 
Run 35 stress 0.1814036 
... Procrustes: rmse 0.03959269  max resid 0.1481975 
Run 36 stress 0.1822264 
Run 37 stress 0.1825052 
Run 38 stress 0.1818223 
Run 39 stress 0.1810379 
... Procrustes: rmse 0.01367068  max resid 0.08576246 
Run 40 stress 0.1837642 
Run 41 stress 0.1823669 
Run 42 stress 0.1815993 
Run 43 stress 0.1863375 
Run 44 stress 0.1810309 
... Procrustes: rmse 0.01329809  max resid 0.08407502 
Run 45 stress 0.1824075 
Run 46 stress 0.1970325 
Run 47 stress 0.1810496 
... Procrustes: rmse 0.02723117  max resid 0.1117951 
Run 48 stress 0.1815987 
Run 49 stress 0.1828642 
Run 50 stress 0.1831933 
Run 51 stress 0.1814147 
... Procrustes: rmse 0.02883545  max resid 0.1054407 
Run 52 stress 0.1830535 
Run 53 stress 0.1812633 
... Procrustes: rmse 0.01538529  max resid 0.08064122 
Run 54 stress 0.1812549 
... Procrustes: rmse 0.01637034  max resid 0.08068132 
Run 55 stress 0.1821578 
Run 56 stress 0.1881487 
Run 57 stress 0.1814143 
... Procrustes: rmse 0.02873776  max resid 0.1051949 
Run 58 stress 0.1873034 
Run 59 stress 0.1810298 
... Procrustes: rmse 0.01306299  max resid 0.08278488 
Run 60 stress 0.1814024 
... Procrustes: rmse 0.03943623  max resid 0.1473544 
Run 61 stress 0.1809421 
... Procrustes: rmse 0.0198368  max resid 0.0992024 
Run 62 stress 0.1810715 
... Procrustes: rmse 0.03980534  max resid 0.1467033 
Run 63 stress 0.182157 
Run 64 stress 0.1833543 
Run 65 stress 0.1810316 
... Procrustes: rmse 0.01158711  max resid 0.07286033 
Run 66 stress 0.1814024 
... Procrustes: rmse 0.03969086  max resid 0.1472131 
Run 67 stress 0.1821248 
Run 68 stress 0.1823823 
Run 69 stress 0.1819351 
Run 70 stress 0.1816261 
Run 71 stress 0.182226 
Run 72 stress 0.1814204 
... Procrustes: rmse 0.03583938  max resid 0.1412563 
Run 73 stress 0.1819335 
Run 74 stress 0.1822758 
Run 75 stress 0.1819161 
Run 76 stress 0.1816036 
Run 77 stress 0.1814157 
... Procrustes: rmse 0.02745011  max resid 0.1035174 
Run 78 stress 0.1810491 
... Procrustes: rmse 0.02698331  max resid 0.1109183 
Run 79 stress 0.1812547 
... Procrustes: rmse 0.01650078  max resid 0.08064451 
Run 80 stress 0.1821955 
Run 81 stress 0.1813262 
... Procrustes: rmse 0.02223434  max resid 0.09864611 
Run 82 stress 0.1808178 
... New best solution
... Procrustes: rmse 0.009380364  max resid 0.04964183 
Run 83 stress 0.1873038 
Run 84 stress 0.1822261 
Run 85 stress 0.182276 
Run 86 stress 0.1814022 
Run 87 stress 0.1832011 
Run 88 stress 0.1881501 
Run 89 stress 0.1907288 
Run 90 stress 0.1838924 
Run 91 stress 0.1818229 
Run 92 stress 0.1808181 
... Procrustes: rmse 0.0002255994  max resid 0.001199249 
... Similar to previous best
*** Best solution repeated 1 times

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)  + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot

library(cowplot)
my_legend <- get_legend(T1_rare_plot)
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.pdf")

Saving 7.29 x 4.51 in image

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.eps")

Saving 7.29 x 4.51 in image

CAP ordination plots rarefied - not used

t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1_cap.pdf")


t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2_cap.pdf")


#G166SG identified as outlier based on plots with it included. Removed to create plot. 
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3_cap.pdf")
ggplot_build(T3_rare$plot)$data

Pairwise adonis testing

ps_pairwiseadonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}

ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)

ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)

ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)


ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)

Pairwise betadispr by treatment, time and mode

ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
                mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
        ## Perform test
                print(anova(mod))
        ## Permutation test for F
                pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
                print(pmod)
                print(boxplot(mod))
}

permute test of disperson

ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value    Pr(>F)    
Groups     2 0.010804 0.0054021  11.752 5.966e-05 ***
Residuals 53 0.024363 0.0004597                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm   Pr(>F)    
Groups     2 0.010804 0.0054021 11.752   1000 0.000999 ***
Residuals 53 0.024363 0.0004597                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Chemical       Hand Non-Treated
Chemical               0.00099900      0.0030
Hand        0.00014786                 0.2198
Non-Treated 0.00224710 0.23207124            
$stats
          [,1]      [,2]      [,3]
[1,] 0.3804321 0.3504290 0.3665005
[2,] 0.3977447 0.3647856 0.3788662
[3,] 0.4169569 0.3742943 0.3945778
[4,] 0.4264042 0.4009722 0.4033251
[5,] 0.4454283 0.4306733 0.4230861

$n
[1] 34 10 12

$conf
          [,1]      [,2]      [,3]
[1,] 0.4091911 0.3562140 0.3834219
[2,] 0.4247227 0.3923745 0.4057337

$out
   G046SG 
0.4892527 

$group
[1] 1

$names
[1] "Chemical"    "Hand"        "Non-Treated"

ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups     2 0.005564 0.0027818  2.6163 0.08286 .
Residuals 51 0.054227 0.0010633                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm  Pr(>F)  
Groups     2 0.005564 0.0027818 2.6163   1000 0.07892 .
Residuals 51 0.054227 0.0010633                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Chemical     Hand Non-Treated
Chemical             0.563437      0.0180
Hand        0.544593               0.2138
Non-Treated 0.016124 0.196261            
$stats
          [,1]      [,2]      [,3]
[1,] 0.3663187 0.3540293 0.3405369
[2,] 0.3828402 0.3788455 0.3696567
[3,] 0.3940453 0.3912101 0.3843606
[4,] 0.4182627 0.4028680 0.3904450
[5,] 0.4630042 0.4051548 0.4046331

$n
[1] 32 11 11

$conf
          [,1]      [,2]      [,3]
[1,] 0.3841516 0.3797660 0.3744573
[2,] 0.4039391 0.4026541 0.3942639

$out
   G070SG    G065SG 
0.5218521 0.5170988 

$group
[1] 1 2

$names
[1] "Chemical"    "Hand"        "Non-Treated"

ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.000037 0.00001841  0.0189 0.9813
Residuals 51 0.049706 0.00097462               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000037 0.00001841 0.0189   1000  0.982
Residuals 51 0.049706 0.00097462                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Chemical    Hand Non-Treated
Chemical             0.95604      0.8262
Hand         0.95042              0.9421
Non-Treated  0.81436 0.93161            
$stats
          [,1]      [,2]      [,3]
[1,] 0.3660073 0.3529787 0.3641228
[2,] 0.3775720 0.3619575 0.3888203
[3,] 0.3938294 0.3833167 0.3943978
[4,] 0.4094830 0.4127647 0.4146763
[5,] 0.4548827 0.4517078 0.4523815

$n
[1] 33 11 10

$conf
          [,1]      [,2]      [,3]
[1,] 0.3850526 0.3591128 0.3814792
[2,] 0.4026063 0.4075206 0.4073165

$out
   G123SG    G125SG 
0.4884281 0.5089781 

$group
[1] 1 2

$names
[1] "Chemical"    "Hand"        "Non-Treated"

ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value  Pr(>F)  
Groups     2 0.005636 0.00281822  3.8869 0.02381 *
Residuals 96 0.069606 0.00072506                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm  Pr(>F)  
Groups     2 0.005636 0.00281822 3.8869   1000 0.01798 *
Residuals 96 0.069606 0.00072506                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          T1        T2     T3
T1           0.2007992 0.0030
T2 0.1962534           0.2238
T3 0.0023253 0.1971546       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3804400 0.3663109 0.3660212
[2,] 0.3977547 0.3828261 0.3776147
[3,] 0.4169648 0.3940558 0.3939165
[4,] 0.4264004 0.4182577 0.4094472
[5,] 0.4454327 0.4630282 0.4548610

$n
[1] 34 32 33

$conf
          [,1]      [,2]      [,3]
[1,] 0.4092027 0.3841594 0.3851612
[2,] 0.4247268 0.4039521 0.4026718

$out
   G046SG    G070SG    G123SG 
0.4892518 0.5218707 0.4883976 

$group
[1] 1 2 3

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df    Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.0019364 0.00096822  2.4191 0.1062
Residuals 30 0.0120071 0.00040024               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.0019364 0.00096822 2.4191   1000 0.1049
Residuals 30 0.0120071 0.00040024                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         T1       T2     T3
T1          0.088911 0.5215
T2 0.099548          0.0579
T3 0.496758 0.059679       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3665168 0.3405704 0.3641416
[2,] 0.3788611 0.3696647 0.3888093
[3,] 0.3945750 0.3843826 0.3943811
[4,] 0.4033384 0.3904236 0.4146615
[5,] 0.4230848 0.4046057 0.4524094

$n
[1] 12 11 10

$conf
          [,1]      [,2]      [,3]
[1,] 0.3834107 0.3744933 0.3814644
[2,] 0.4057392 0.3942719 0.4072979

$out
numeric(0)

$group
numeric(0)

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.001904 0.00095187  0.5916   0.56
Residuals 29 0.046661 0.00160901               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.001904 0.00095187 0.5916   1000 0.5734
Residuals 29 0.046661 0.00160901                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        T1      T2     T3
T1         0.31868 0.3896
T2 0.27871         0.9451
T3 0.35755 0.94959       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3504614 0.3540454 0.3530162
[2,] 0.3647810 0.3788398 0.3619277
[3,] 0.3742906 0.3912157 0.3833179
[4,] 0.4009630 0.4028635 0.4127718
[5,] 0.4306622 0.4051523 0.4516826

$n
[1] 10 11 11

$conf
          [,1]      [,2]      [,3]
[1,] 0.3562127 0.3797711 0.3590964
[2,] 0.3923686 0.4026603 0.4075394

$out
   G065SG    G125SG 
0.5171023 0.5089702 

$group
[1] 2 3

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Herbicide")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     4 0.002840 0.00070997   1.074 0.3791
Residuals 51 0.033714 0.00066106               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq     F N.Perm Pr(>F)
Groups     4 0.002840 0.00070997 1.074   1000 0.3886
Residuals 51 0.033714 0.00066106                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                    Atrazine-Mesotrione  Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione                     0.925075   0.966034   0.086913      0.3646
Dicamba                        0.930200            0.978022   0.110889      0.3996
Glyphosate                     0.967636 0.975250              0.188811      0.5365
Handweeded                     0.079155 0.109142   0.175310                 0.2298
Non-Treated                    0.369597 0.406696   0.506436   0.232071            
$stats
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3667119 0.3551354 0.3573486 0.3504290 0.3665005
[2,] 0.3836426 0.3774236 0.3790031 0.3647856 0.3788662
[3,] 0.4069157 0.4021368 0.3897409 0.3742943 0.3945778
[4,] 0.4088518 0.4231790 0.4188086 0.4009722 0.4033251
[5,] 0.4393221 0.4426598 0.4302367 0.4306733 0.4230861

$n
[1] 11 12 11 10 12

$conf
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3949063 0.3812674 0.3707780 0.3562140 0.3834219
[2,] 0.4189251 0.4230061 0.4087038 0.3923745 0.4057337

$out
   G046SG 
0.4804066 

$group
[1] 3

$names
[1] "Atrazine-Mesotrione" "Dicamba"             "Glyphosate"          "Handweeded"          "Non-Treated"        

ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Herbicide")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     4 0.002843 0.00071079  0.6541 0.6268
Residuals 49 0.053249 0.00108672               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     4 0.002843 0.00071079 0.6541   1000 0.6404
Residuals 49 0.053249 0.00108672                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                    Atrazine-Mesotrione  Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione                     0.742258   0.675325   0.735265      0.4426
Dicamba                        0.727584            0.313686   0.971029      0.0739
Glyphosate                     0.640120 0.282922              0.399600      0.6204
Handweeded                     0.734970 0.963589   0.382841                 0.1918
Non-Treated                    0.393767 0.097876   0.590421   0.196261            
$stats
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3474546 0.3612943 0.3592504 0.3540293 0.3405369
[2,] 0.3649680 0.3733844 0.3708470 0.3788455 0.3696567
[3,] 0.3810085 0.3939371 0.3751855 0.3912101 0.3843606
[4,] 0.4036468 0.4158683 0.4004167 0.4028680 0.3904450
[5,] 0.4120948 0.4539607 0.4344526 0.4051548 0.4046331

$n
[1] 11 11 10 11 11

$conf
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3625824 0.3736982 0.3604113 0.3797660 0.3744573
[2,] 0.3994346 0.4141759 0.3899598 0.4026541 0.3942639

$out
   G070SG    G065SG 
0.5064339 0.5170988 

$group
[1] 1 4

$names
[1] "Atrazine-Mesotrione" "Dicamba"             "Glyphosate"          "Handweeded"          "Non-Treated"        

ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     4 0.003957 0.00098924  1.0438 0.3944
Residuals 49 0.046441 0.00094777               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     4 0.003957 0.00098924 1.0438   1000 0.3806
Residuals 49 0.046441 0.00094777                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                    Atrazine-Mesotrione  Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione                     0.369630   0.710290   0.311688      0.1049
Dicamba                        0.355725            0.193806   0.764236      0.5584
Glyphosate                     0.696633 0.191707              0.225774      0.0330
Handweeded                     0.305819 0.741563   0.222883                 0.9351
Non-Treated                    0.107926 0.561747   0.033811   0.931610            
$stats
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3541179 0.3627381 0.3520382 0.3529787 0.3641228
[2,] 0.3639605 0.3754509 0.3633215 0.3619575 0.3888203
[3,] 0.3735777 0.3869820 0.3798037 0.3833167 0.3943978
[4,] 0.3928099 0.3954294 0.3890832 0.4127647 0.4146763
[5,] 0.4082574 0.4037548 0.4035468 0.4517078 0.4523815

$n
[1] 12 11 10 11 10

$conf
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3604193 0.3774645 0.3669321 0.3591128 0.3814792
[2,] 0.3867361 0.3964996 0.3926753 0.4075206 0.4073165

$out
   G131SG    G123SG    G125SG 
0.4392663 0.4725323 0.5089781 

$group
[1] 1 2 4

$names
[1] "Atrazine-Mesotrione" "Dicamba"             "Glyphosate"          "Handweeded"          "Non-Treated"        

ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df    Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.0028576 0.00142878  2.0791 0.1439
Residuals 28 0.0192418 0.00068721               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.0028576 0.00142878 2.0791   1000 0.1678
Residuals 28 0.0192418 0.00068721                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         T1       T2     T3
T1          0.288711 0.0709
T2 0.260682          0.4146
T3 0.071566 0.402183       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3573164 0.3592595 0.3520475
[2,] 0.3789992 0.3708523 0.3633372
[3,] 0.3897552 0.3751867 0.3798028
[4,] 0.4188308 0.4004023 0.3890938
[5,] 0.4302344 0.4344627 0.4035369

$n
[1] 11 10 10

$conf
          [,1]      [,2]      [,3]
[1,] 0.3707798 0.3604224 0.3669337
[2,] 0.4087305 0.3899511 0.3926719

$out
   G046SG 
0.4803861 

$group
[1] 1

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df    Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.0020678 0.00103388  1.0923  0.348
Residuals 31 0.0293417 0.00094651               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.0020678 0.00103388 1.0923   1000 0.3716
Residuals 31 0.0293417 0.00094651                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         T1       T2     T3
T1          0.658342 0.0500
T2 0.613621          0.4595
T3 0.059278 0.439855       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3667106 0.3474538 0.3540744
[2,] 0.3836487 0.3649640 0.3639284
[3,] 0.4069173 0.3810218 0.3735603
[4,] 0.4088534 0.4036548 0.3928858
[5,] 0.4393202 0.4120929 0.4082313

$n
[1] 11 11 12

$conf
          [,1]      [,2]      [,3]
[1,] 0.3949101 0.3625900 0.3603526
[2,] 0.4189245 0.3994537 0.3867680

$out
   G070SG    G131SG 
0.5064274 0.4391984 

$group
[1] 2 3

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df    Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.0004756 0.00023778  0.2821 0.7561
Residuals 31 0.0261290 0.00084287               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.0004756 0.00023778 0.2821   1000 0.7572
Residuals 31 0.0261290 0.00084287                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        T1      T2     T3
T1         0.82118 0.4815
T2 0.82200         0.6314
T3 0.46827 0.62758       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3551246 0.3612819 0.3627463
[2,] 0.3774302 0.3733980 0.3754452
[3,] 0.4021407 0.3938902 0.3869939
[4,] 0.4231764 0.4158655 0.3954270
[5,] 0.4426504 0.4539101 0.4037696

$n
[1] 12 11 11

$conf
          [,1]      [,2]      [,3]
[1,] 0.3812756 0.3736592 0.3774747
[2,] 0.4230059 0.4141212 0.3965130

$out
   G123SG 
0.4725292 

$group
[1] 3

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.001904 0.00095187  0.5916   0.56
Residuals 29 0.046661 0.00160901               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     2 0.001904 0.00095187 0.5916   1000 0.5774
Residuals 29 0.046661 0.00160901                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        T1      T2     T3
T1         0.31668 0.3566
T2 0.27871         0.9431
T3 0.35755 0.94959       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3504614 0.3540454 0.3530162
[2,] 0.3647810 0.3788398 0.3619277
[3,] 0.3742906 0.3912157 0.3833179
[4,] 0.4009630 0.4028635 0.4127718
[5,] 0.4306622 0.4051523 0.4516826

$n
[1] 10 11 11

$conf
          [,1]      [,2]      [,3]
[1,] 0.3562127 0.3797711 0.3590964
[2,] 0.3923686 0.4026603 0.4075394

$out
   G065SG    G125SG 
0.5171023 0.5089702 

$group
[1] 2 3

$names
[1] "T1" "T2" "T3"

ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")

Analysis of Variance Table

Response: Distances
          Df    Sum Sq    Mean Sq F value Pr(>F)
Groups     2 0.0019364 0.00096822  2.4191 0.1062
Residuals 30 0.0120071 0.00040024               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.0019364 0.00096822 2.4191   1000 0.0989 .
Residuals 30 0.0120071 0.00040024                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         T1       T2     T3
T1          0.090909 0.5145
T2 0.099548          0.0639
T3 0.496758 0.059679       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3665168 0.3405704 0.3641416
[2,] 0.3788611 0.3696647 0.3888093
[3,] 0.3945750 0.3843826 0.3943811
[4,] 0.4033384 0.3904236 0.4146615
[5,] 0.4230848 0.4046057 0.4524094

$n
[1] 12 11 10

$conf
          [,1]      [,2]      [,3]
[1,] 0.3834107 0.3744933 0.3814644
[2,] 0.4057392 0.3942719 0.4072979

$out
numeric(0)

$group
numeric(0)

$names
[1] "T1" "T2" "T3"

ps_betadispr(ps_rare_sub, groupingvar = "Herbicide")

Analysis of Variance Table

Response: Distances
           Df   Sum Sq    Mean Sq F value Pr(>F)
Groups      4 0.003022 0.00075548  0.8548 0.4926
Residuals 159 0.140523 0.00088379               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups      4 0.003022 0.00075548 0.8548   1000 0.5185
Residuals 159 0.140523 0.00088379                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                    Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione                     0.29870    0.73926    0.46454      0.8711
Dicamba                         0.27184            0.11988    0.90010      0.1299
Glyphosate                      0.73022 0.10681               0.28272      0.7942
Handweeded                      0.45118 0.88577    0.27400                 0.3097
Non-Treated                     0.87491 0.12162    0.79878    0.31930            
$stats
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3635263 0.3726987 0.3644354 0.3712730 0.3653358
[2,] 0.3806232 0.3917390 0.3838596 0.3907178 0.3894408
[3,] 0.3980967 0.4079200 0.3980853 0.3985758 0.4017338
[4,] 0.4199957 0.4233531 0.4072361 0.4212779 0.4126745
[5,] 0.4700021 0.4662951 0.4375529 0.4411483 0.4390542

$n
[1] 34 34 31 32 33

$conf
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3874280 0.3993536 0.3914516 0.3900402 0.3953436
[2,] 0.4087653 0.4164864 0.4047190 0.4071115 0.4081241

$out
   G070SG    G123SG    G046SG    G076SG    G065SG    G125SG    G151SG 
0.5184955 0.4995460 0.4894611 0.4554918 0.5311635 0.5552741 0.4631531 

$group
[1] 1 2 3 3 4 4 5

$names
[1] "Atrazine-Mesotrione" "Dicamba"             "Glyphosate"          "Handweeded"          "Non-Treated"        

ps_betadispr(ps_rare_sub, groupingvar = "Mode")

Analysis of Variance Table

Response: Distances
           Df   Sum Sq    Mean Sq F value Pr(>F)
Groups      2 0.001682 0.00084096  0.9478 0.3897
Residuals 161 0.142847 0.00088725               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups      2 0.001682 0.00084096 0.9478   1000 0.3746
Residuals 161 0.142847 0.00088725                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Chemical    Hand Non-Treated
Chemical             0.97502      0.1279
Hand         0.97790              0.2827
Non-Treated  0.13961 0.31930            
$stats
          [,1]      [,2]      [,3]
[1,] 0.3678724 0.3712730 0.3653358
[2,] 0.3914997 0.3907178 0.3894408
[3,] 0.4069683 0.3985758 0.4017338
[4,] 0.4198506 0.4212779 0.4126745
[5,] 0.4591917 0.4411483 0.4390542

$n
[1] 99 32 33

$conf
          [,1]      [,2]      [,3]
[1,] 0.4024663 0.3900402 0.3953436
[2,] 0.4114703 0.4071115 0.4081241

$out
   G046SG    G070SG    G087SG    G123SG    G131SG    G065SG    G125SG    G151SG 
0.4911052 0.5250529 0.4726503 0.4968193 0.4692099 0.5311635 0.5552741 0.4631531 

$group
[1] 1 1 1 1 1 2 2 3

$names
[1] "Chemical"    "Hand"        "Non-Treated"

ps_betadispr(ps_rare_sub, groupingvar = "Time")

Analysis of Variance Table

Response: Distances
           Df   Sum Sq    Mean Sq F value Pr(>F)
Groups      2 0.001119 0.00055970  0.6361 0.5307
Residuals 161 0.141654 0.00087984               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups      2 0.001119 0.00055970 0.6361   1000 0.5455
Residuals 161 0.141654 0.00087984                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        T1      T2     T3
T1         0.29271 0.3367
T2 0.29177         0.9800
T3 0.31167 0.98337       
$stats
          [,1]      [,2]      [,3]
[1,] 0.3751393 0.3651628 0.3656380
[2,] 0.3983422 0.3889047 0.3843976
[3,] 0.4102601 0.3985027 0.4009712
[4,] 0.4266630 0.4186156 0.4162396
[5,] 0.4483672 0.4529592 0.4619116

$n
[1] 56 54 54

$conf
          [,1]      [,2]      [,3]
[1,] 0.4042806 0.3921145 0.3941248
[2,] 0.4162397 0.4048909 0.4078175

$out
   G046SG    G065SG    G070SG    G087SG    G123SG    G125SG 
0.4859161 0.5363026 0.5212376 0.4661354 0.4943452 0.5607126 

$group
[1] 1 2 2 2 3 3

$names
[1] "T1" "T2" "T3"

beta_boxplot<-function (physeq, method = "bray", group) 
{
  require("phyloseq")
  require("ggplot2")
  group2samp <- list()
  group_list <- get_variable(sample_data(physeq), group)
  for (groups in levels(group_list)) {
    target_group <- which(group_list == groups)
    group2samp[[groups]] <- sample_names(physeq)[target_group]
  }
  beta_div_dist <- phyloseq::distance(physeq = physeq, method = method)
  beta_div_dist <- as(beta_div_dist, "matrix")
  dist_df <- data.frame()
  counter <- 1
  for (groups in names(group2samp)) {
    sub_dist <- beta_div_dist[group2samp[[groups]], group2samp[[groups]]]
    no_samp_col <- ncol(sub_dist)
    no_samp_row <- nrow(sub_dist)
    for (cols in seq(no_samp_col)) {
      if (cols > 1) {
        for (rows in seq((cols - 1))) {
          dist_df[counter, "sample_pair"] <- paste0(colnames(sub_dist)[cols], 
            "-", rownames(sub_dist)[rows])
          dist_df[counter, "group"] <- groups
          dist_df[counter, "beta_div_method"] <- method
          dist_df[counter, "beta_div_value"] <- sub_dist[rows, 
            cols]
          counter = counter + 1
        }
      }
    }
  }
  plot_boxplot <- ggplot(data = dist_df, aes(x = group, y = beta_div_value, 
    color = group)) + geom_boxplot(outlier.shape = NA) + 
    geom_jitter() + theme_bw() + xlab(group) + ylab(method) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, 
      hjust = 1))
  list_Out <- list(data = dist_df, plot = plot_boxplot)
  return(list_Out)
}

box and whisker plots of pairwise distance within group distances

#remotes::install_github("antonioggsousa/micrUBIfuns")
#library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T1_beta_plot

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

my_legend <- get_legend(T1_beta_plot)

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.eps")

Saving 7.29 x 4.51 in image

T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none") 
T1_beta_plot

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.02391 0.005978     5.4 0.000334 ***
Residuals   282 0.31223 0.001107                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T1_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)

$group
                                        diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione     -0.001507071 -0.018186144  0.015172002 0.9991606
Glyphosate-Atrazine-Mesotrione   0.001896970 -0.015523754  0.019317693 0.9982525
Handweeded-Atrazine-Mesotrione  -0.023576431 -0.041939486 -0.005213376 0.0044684
Non-Treated-Atrazine-Mesotrione -0.013193939 -0.029873012  0.003485134 0.1934686
Glyphosate-Dicamba               0.003404040 -0.013275033  0.020083113 0.9805780
Handweeded-Dicamba              -0.022069360 -0.039730381 -0.004408340 0.0061800
Non-Treated-Dicamba             -0.011686869 -0.027589741  0.004216003 0.2601887
Handweeded-Glyphosate           -0.025473401 -0.043836456 -0.007110346 0.0015994
Non-Treated-Glyphosate          -0.015090909 -0.031769982  0.001588164 0.0971493
Non-Treated-Handweeded           0.010382492 -0.007278529  0.028043512 0.4896459

T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T2_beta_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)

             Df Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.0328 0.008212   5.471 0.000303 ***
Residuals   260 0.3902 0.001501                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T2_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)

$group
                                        diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione      0.004587879 -0.015705828  0.024881585 0.9716343
Glyphosate-Atrazine-Mesotrione  -0.011421549 -0.032812994  0.009969896 0.5850476
Handweeded-Atrazine-Mesotrione   0.007709091 -0.012584615  0.028002797 0.8348219
Non-Treated-Atrazine-Mesotrione -0.022036364 -0.042330070 -0.001742657 0.0257655
Glyphosate-Dicamba              -0.016009428 -0.037400872  0.005382017 0.2426725
Handweeded-Dicamba               0.003121212 -0.017172494  0.023414919 0.9933115
Non-Treated-Dicamba             -0.026624242 -0.046917949 -0.006330536 0.0034255
Handweeded-Glyphosate            0.019130640 -0.002260805  0.040522085 0.1039333
Non-Treated-Glyphosate          -0.010614815 -0.032006260  0.010776630 0.6518171
Non-Treated-Handweeded          -0.029745455 -0.050039161 -0.009451748 0.0007049

T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide") 
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500",  "#00BF7D", "#00B0F6", "#E76BF3")) + 
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)

             Df Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.0477 0.011924   10.39 7.93e-08 ***
Residuals   261 0.2994 0.001147                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T3_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)

$group
                                         diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione      0.0183969697  0.001410694  0.035383246 0.0263568
Glyphosate-Atrazine-Mesotrione  -0.0007666667 -0.018752976  0.017219643 0.9999574
Handweeded-Atrazine-Mesotrione   0.0276454545  0.010659179  0.044631731 0.0001130
Non-Treated-Atrazine-Mesotrione  0.0315000000  0.013513690  0.049486310 0.0000250
Glyphosate-Dicamba              -0.0191636364 -0.037864911 -0.000462362 0.0415656
Handweeded-Dicamba               0.0092484848 -0.008493102  0.026990072 0.6075977
Non-Treated-Dicamba              0.0131030303 -0.005598244  0.031804305 0.3068486
Handweeded-Glyphosate            0.0284121212  0.009710847  0.047113396 0.0003917
Non-Treated-Glyphosate           0.0322666667  0.012652605  0.051880729 0.0000917
Non-Treated-Handweeded           0.0038545455 -0.014846729  0.022555820 0.9798083

library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 5, height = 10)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.eps", width = 5, height = 10)

Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.

T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"


beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)

beta_anova<-function(data, Herbicide = "Herbicide"){
  df_sub<- data %>% filter(group == Herbicide)
  mod<-aov(beta_div_value ~ Time, data = df_sub)
  print(summary(mod))
  print(TukeyHSD(mod, "Time"))
  boxplot(df_sub$beta_div_value ~ df_sub$Time)
}

beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")

             Df  Sum Sq Mean Sq F value   Pr(>F)    
Time          2 0.02547 0.01274   17.94 9.08e-08 ***
Residuals   163 0.11573 0.00071                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff          lwr          upr     p adj
T2-T1 -0.01621212 -0.027718988 -0.004705255 0.0030431
T3-T1  0.01578788  0.003603568  0.027972190 0.0071575
T3-T2  0.03200000  0.019331357  0.044668643 0.0000000

beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")

             Df  Sum Sq  Mean Sq F value Pr(>F)  
Time          2 0.01713 0.008567   4.195 0.0168 *
Residuals   152 0.31041 0.002042                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff           lwr        upr     p adj
T2-T1  0.02391582  0.0024158024 0.04541585 0.0251783
T3-T1  0.02231582  0.0008158024 0.04381585 0.0399463
T3-T2 -0.00160000 -0.0219967123 0.01879671 0.9811772

beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")

             Df  Sum Sq  Mean Sq F value Pr(>F)
Time          2 0.00273 0.001366   0.977  0.378
Residuals   173 0.24171 0.001397               
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
              diff         lwr         upr     p adj
T2-T1 -0.001274747 -0.01740797 0.014858477 0.9809505
T3-T1 -0.009002020 -0.02513524 0.007131204 0.3864971
T3-T2 -0.007727273 -0.02457788 0.009123330 0.5252045

beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
Time          2 0.02773 0.013867   11.41 2.22e-05 ***
Residuals   173 0.21026 0.001215                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
              diff         lwr          upr     p adj
T2-T1 -0.007369697 -0.02308574  0.008346347 0.5100816
T3-T1 -0.028906061 -0.04395303 -0.013859094 0.0000310
T3-T2 -0.021536364 -0.03658333 -0.006489397 0.0025367

beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
Time          2 0.02597 0.012985    14.9 1.33e-06 ***
Residuals   142 0.12374 0.000871                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff         lwr          upr     p adj
T2-T1 -0.02068822 -0.03474249 -0.006633939 0.0018742
T3-T1 -0.03156970 -0.04562397 -0.017515421 0.0000012
T3-T2 -0.01088148 -0.02562173  0.003858768 0.1909546

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 


T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot 
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T1_beta_chemical_combined_plot

Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

my_legend <- get_legend(T1_beta_chemical_combined_plot)

Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.eps")

Saving 7.29 x 4.51 in image

T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot

Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot 
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T2_beta_chemical_combined_plot

T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot 
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T3_beta_chemical_combined_plot

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)

Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.eps", width = 5, height = 10)

T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data

T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"

m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)


m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)

m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)


beta_div_chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)

beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  #wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta_multicomparison.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta_multicomparison.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta_multicomparison.pdf")

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta_multi_comparison.pdf", width = 7, height = 10)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Taxon_barplot_t1.pdf")

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t2.pdf")

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t3.pdf")

Combined herbicide and time bar plot for exploration

sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare) > 50, ps_rare)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)

Linear modeling of abundant taxa.

sig_mods

NULL